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Abstract 

Bounds are developed for the condition number of the linear finite element equations of an 
anisotropic diffusion problem with arbitrary meshes. They depend on three factors. The first, 
factor proportional to a power of the number of mesh elements, represents the condition number 
of the linear finite element equations for the Laplacian operator on a uniform mesh. The other 
two factors arise from the mesh nonuniformity viewed in the Euclidean metric and in the metric 
defined by the diffusion matrix. The new bounds reveal that the conditioning of the finite 
element equations with adaptive anisotropic meshes is much better than what is commonly 
feared. Diagonal scaling for the linear system and its effects on the conditioning are also studied. 
It is shown that the Jacobi preconditioning, which is an optimal diagonal scaling for a symmetric 
positive definite sparse matrix, can eliminate the effects of mesh nonuniformity viewed in the 
Euclidean metric and reduce those effects of the mesh viewed in the metric defined by the 
diffusion matrix. Tight bounds on the extreme eigenvalues of the stiffness and mass matrices are 
obtained. Numerical examples are given. 
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1 Introduction 



It has been amply demonstrated that significant improvements in accuracy can be gained when an 
appropriately chosen anisotropic mesh is used for the numerical solution of problems exhibiting 
anisotropic features. However, there exists a general concern in the scientific computing community 
that an anisotropic mesh, which can contain elements of large aspect ratio and small volume, 
may lead to ill-conditioned linear systems and this could outweigh the accuracy and efficiency 
improvements gained by anisotropic mesh adaptation. For isotropic mesh adaptation, Bank and 
Scott [2] (also see Brenner and Scott [3]) show that after proper diagonal scaling, the condition 
number of finite element equations with an adaptive mesh is essentially the same as that for a 
uniform mesh. Unfortunately, this result does not apply to anisotropic meshes nor to problems with 
anisotropic diffusion. 

For problems with anisotropic diffusion and arbitrary meshes, several estimates have been 
developed for the extreme eigenvalues of the stiffness matrix. For example, Fried [7] shows that the 
largest eigenvalue of the stiffness matrix is bounded by the largest eigenvalues of element stiffness 
matrices. Shewchuk [18 obtains sharp bounds on largest eigenvalues of element stiffness matrices 
for linear triangular and tetrahedral finite elements. More recently, Du et al. [5] develop a bound 
that can be viewed as a generalization of Shewchuk's result to general dimensions and simplicial 
finite elements. 

Estimation of the smallest eigenvalue for the general case appears to be more challenging. 
Standard estimates (e.g., see Ern and Guermond |6j) are linearly proportional to the volume of 
the smallest mesh element, which is typically too pessimistic for nonuniform meshes. Moreover, 
Apel JTJ Sect. 4.3.3] shows that the order of the smallest eigenvalue of the stiffness matrix for a 
specific, specially designed anisotropic mesh is the same as for a uniform mesh. As a matter of fact, 
coefficient adaptive anisotropic meshes can even improve the conditioning for partial differential 
equations (PDEs) with anisotropic diffusion coefficients, as observed by D'Azavedo et al. [I] and 
Shewchuk |18[ Sect. 3.2]. A noticeable approach for obtaining sharper bounds for the smallest 
eigenvalue is proposed by Fried |Z]. The approach employs a continuous generalized eigenvalue 
problem with an auxiliary density function and its key is to find a lower bound for the smallest 
eigenvalue of the continuous problem. Bounds for the smallest eigenvalue of the stiffness matrix 
obtained with Fried's approach are valid for general meshes in any dimension but in d > 3 dimensions 
they are less sharp than those obtained in this paper. 

The objective of this paper is threefold. First, we develop tight bounds on the extreme eigenvalues 
and the condition number of the stiffness matrix for a general diffusion problem with an arbitrary 
anisotropic mesh. No assumption on the shape or size of mesh elements is made in the development. 
Our upper bound on the largest eigenvalue can be also expressed in terms of mesh nonuniformity 
viewed in the metric tensor defined by the diffusion matrix (which will hereafter be referred to as 
the mesh D -nonuniformity) . It is comparable to those of Shewchuk |18] and Du et al. [S] but is 
expressed as a sum of patch-wise terms instead of element-wise terms as in the aforementioned 
references. The patch-wise nature gives a sharper bound and makes it more convenient to use in the 
development of diagonal scaling preconditioners. To obtain lower bounds for the smallest eigenvalue 
of the stiffness matrix we extend Bank and Scott's result [2] to arbitrary meshes. This generalization 
is not trivial and special effort has to be made to deal with the arbitrariness of the mesh. Along the 
way we establish anisotropic upper and lower bounds on the extreme eigenvalues of the mass matrix 
which are much tighter than estimates available in the literature. 

The second objective of the paper is to provide a clear geometric interpretation for the obtained 
bounds on the condition number of the stiffness matrix. These bounds are shown to depend on three 
factors. The first factor is proportional to a power of the number of mesh elements and represents 
the condition number of the stiffness matrix for the linear finite element approximation of the 
Laplacian operator on a uniform mesh. The other two factors arise from the mesh nonuniformity in 
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volume measured in the Euclidean metric (which will be referred to the mesh volume-nonuniformity) 
and from the mesh D-nonuniformity. 

The third objective is to study diagonal scaling for the finite element linear system and its 
effects on the conditioning. We focus on the scaling with the diagonal entries of the matrix (Jacobi 
preconditioning) since it is an optimal diagonal scaling for a symmetric positive definite sparse 
matrix |1U| Corollary 7.6 and the following]. We show that the Jacobi preconditioning can eliminate 
the effects of the mesh volume-nonuniformity and improve those caused by the mesh D-nonuniformity, 
thus significantly reducing the effects of the mesh irregularity on the conditioning. From the practical 
point of view, this result indicates that a simple diagonal preconditioning can effectively transform 
the stiffness matrix into a matrix which has a comparable condition number as the one with a 
uniform mesh. 

The outline of the paper is as follows. Section[2]briefly describes a linear finite element discretization 
of a general anisotropic diffusion problem. Estimation of the extreme eigenvalues and the condition 
number of the mass matrix is given in Sect. [3} Section [4] deals with the estimation of the largest 
eigenvalue of the stiffness matrix. Bounds on the smallest eigenvalue and the condition number of 
the stiffness matrix and the effects of diagonal scaling are investigated in Sect. [5j A selection of 
examples in Sect. [6] provides a numerical validation for the theoretical findings. Finally, conclusions 
and further remarks are given in Sect. [7} 

2 Linear finite element approximation 

We consider the boundary value problem (BVP) of a general diffusion differential equation in the 
form 



where ft is a simply connected polygonal or polyhedral domain in M. d (d > 1) and D = D(a:) is the 
diffusion matrix. We assume that D is symmetric and positive definite and there exist two positive 
constants d m ; n and d max such that 



where the less-than-or-equal sign means that the difference between the right-hand side and left-hand 
side terms is positive semidefinite. 

We are interested in the linear finite element solution of BVP ([I]) . Assume that an affine family 
{7~h} of simplicial decompositions of ft is given and denote the associated linear finite element space 
by V h C Ho(Q). A linear finite element solution Uh £ V h to BVP Q is defined by 




V • (DVn) = /, in ft, 
= 0, on <9ft 



(1) 



dminl < D(cc) < (f max I, \/x G ft, 



(2) 




or 




Since both Vu^ and Vvh are constant on K, we can rewrite the above equation as 




(3) 



where Dk is the integral average of D over K, i.e. 




(4) 
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In practice, the integrals in ^ and Q have to be approximated numerically via a quadrature rule. 
Although this will change the definition of D_r- and the right-hand side term of Q slightly, the 
procedure and the results in this paper will remain valid for this situation. 

Finite element equation ^ can be expressed in a matrix form. Denoting the numbers of elements 
and interior vertices of Th by N and N v % and assuming that the vertices are ordered in such a way 
that the first N v i vertices are the interior vertices, we have 

V h = span{0i, . . .,4> Nvi }, 

3 

where <pj is the linear basis function associated with the j th vertex. In ^ and hereafter, we use the 
sum J2j with the index j ranging over all interior vertices, i.e., 

iv ■ 

E = E 

Substituting ([5j into (|3j and taking = (pi (i = 1, . . . , N v i), we obtain the linear algebraic system 

Au = f, 

where u = [m, ■ ■ . , UN vi ] T an d the stiffness matrix A and the right-hand side term / are given by 

A ij= i,j = l,...,N vi , (6) 

KeT h 

fi = ^ f<t>i dx > i = l,..., N vi , 

and V0i|x and V<j)j\K denote the restriction of V</>j and V0j on K. Our main goal is to estimate 
the condition number of stiffness matrix A. 



3 Mass matrix 

To start with, we consider the element mass matrix B for the reference element K (which is assumed 
to be unitary, i.e., \K\ = 1), 

B = (bij), bi tj = J $i4>j d£, i, j = 1, . . . ,d+ 1, 

where <jVs are the linear basis functions associated with the vertices of K. Matrix B is symmetric 
and positive definite. Moreover, for the standard linear finite elements considered in this paper we 
have 

1 KB<J-I. 



(d + l)(d + 2) - ~ d+1 

3.1 Condition number of the mass matrix 

Consider the global mass matrix 



B = (Bi j), B id = / 4>i(f)j dx, i,j = 1 
Jn 



! ■ • • ) N v i . 
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Notice that 

B " - L * dx - x L *> dx = x ( d + i^' + 2) - hMtyy < 7 > 

where coj is the element patch associated with the j th vertex and \ujj\ is its volume. 

The following theorem gives lower and upper bounds on the condition number of the mass matrix 
for any dimension and any mesh. 

Theorem 3.1 (Condition number of the mass matrix). The condition number of the mass matrix 
for the linear finite elements on a simplicial mesh is bounded by 

max; B,-, , , .maxjBjj , . 

3 —^- < k(B) <(d + 2 3 —^-. (8 

min,- Bjj mmj Bjj 

Proof. For an element K, let uk be the restriction of the vector u on K and B% the element mass 
matrix. Then, 



u T 



Bu= UkB k u k = Y \ K I u^Bu K < Y \ k \\\ u k\\1- 

KeT h KET h d+ KeT h 

Rearranging the sum on the right-hand side according to the vertices and using ([T]), 

T 1 1 X — ^ i iii 1 1 2 ^ X — ^ 2 i i ^ X — ^ 2 

u Bu < > \K\ \\uk U = r > Ua uJ = — - — > u,:Bjj, 

"rf+l/i 111 112 d + l^ 3 1 31 2 ^ 3 33 

which implies 

A m ax(£) < ——max Bjj. 



Similarly, we have 



u T Bu > — ^ u2 jBjj and A m j n (i?) > -min Bjj. 



2 ^ -j jj / - 2 

j 



Moreover, it is easy to show that 



Amax(-B) > max Bjj and \ m m{B) < min Bjj. 

j j 



Combining the above estimates gives 

max Bjj < A m a X (-B) < - 1" 2 max J3,-,- , 
j 2 j 

-min Bjj < A min (5) < mni-B^, 

from which estimate ([9| follows. □ 
From ([7]), Theorem 3.1 implies 

U "'* vl < k(B) < (^ + 2)^4, (9) 



where |w m ax| = max.,- \uj\ and |w m in| = mrn il w il- If the mesh is quasi-uniform, this means 
|w ma x|/|w m m| = 0(1) and k(B) = O(l). On the other hand, |w max |/|cj m i n | and k{B) can become 
large for nonuniform meshes. 
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3.2 Relation to the estimates in the literature 



If we denote the maximum number of mesh elements in a patch by p m ax, then 

| -K"min | ^ 1 10 j | 5; Pmax | -^max | , \/j = 1 , . . . , N v { 

and estimate ^ implies 

K (B)<(d + 2) Pm J-§^\, (10) 

I min | 

which is the bound obtained by Fried [7J inequality (24)]. 
Moreover, for an isotropic mesh, 

\K I nc h d \K ■ I oc h d ■ 

l-^maxl ^ "max' K v min| "-min' 

where /i m ax and h m \ n are the largest and smallest element diameters. Substituting this into bound 



(10) gives 

<B)<C(^) . (ID 

fin 



which is precisely the standard estimate found in the literature (e.g., [SJ Rem. 9.10]). 

For anisotropic meshes, on the other hand, the new estimate (|9|) is much tighter than both Fried's 



estimate (10) and the standard estimate (11), since large |if m ax|/|-^mm| and h max /h m i n do not 
necessarily imply large |w m ax|/|w mm | for those meshesQ Furthermore, estimate ^ also provides a 



tight lower bound, which is not available with ( 10 ) and (111 



3.3 Diagonal scaling for the mass matrix 

It is known |10| Corollary 7.6 and the following] that for a symmetric positive definite sparse matrix, 
scaling by its diagonal entries (Jacobi preconditioning) is an optimal diagonal preconditioning (up 
to a constant depending on the maximum number of non-zeroes per column and row of the matrix). 
We are interested in a bound on the condition number after such preconditioning. 



For a diagonal scaling S = (sj), similarly to Theorem 3.1 we obtain 



max,- 2 Ba -. , max, s - 2 Baa 

3 J _ 2 n < KiS^BS- 1 ) < (d + 2) 3 J _ 2 " 

mm i s i B n mm jSj !'>,, 

and, for the Jacobi preconditioning s 2 = Bjj, we have arrived at the following theorem by Wathen jTU] 
who studies the effects of the diagonal scaling on the condition number of the Galerkin mass matrix. 

Theorem 3.2 ([19, Table 1]). The condition number of the Jacobi preconditioned Galerkin mass 
matrix with a simplicial mesh has a mesh-independent bound 

KiS^BS- 1 ) <d + 2. 



Theorems 3.1 and 3.2 show that the mesh volume-nonuniformity has a significant effect on 
the condition number of the mass matrix and this effect is completely eliminated by the Jacobi 
preconditioning. 

As we will see later in Sect. |5j diagonal scaling plays a similar role in reducing the effects of mesh 
nonuniformity on the condition number of the stiffness matrix. 



1 For example, meshes in Fig. [2]have | K m&x \ / \ A' min — > oo but |i^ max |/|^mm| = 0(1). 
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4 Largest eigenvalue of the stiffness matrix 

The following lemma is valid for any dimension. 

Lemma 4.1 (Largest eigenvalue). The largest eigenvalue of the stiffness matrix A = (Aij) for the 
linear finite element approximation of BVP ([IJ is bounded by 

max Ajj < A max (^4) < (d + l)maxA,-,-. (12) 

3 3 

The largest eigenvalue of the diagonally (Jacobi) preconditioned stiffness matrix S^^^AS^ 1 has a 
mesh-independent bound 

1< A max (S'- 1 AS- 1 )<d + l. (13) 
Proof. First, recall that for any symmetric positive semidefmite matrix M, 

v T Mw < - (y T Mv + w T Mw S j , Vt>, w e R d+1 . 

Then, using the local indices on K and the definition of Ajj from ([6| and rearranging the sum 
according to the vertices, we have 



u T Au = / Vuh ■ DVti/j dx 
Jn 



d+l 

= J2 \ K \ E KV^J.Djf^Vfe) 

KeT h iK,3K=l 

d+l 

<(d+i)j2 \ K \ E («kVfe)'%(%v^ K ; 

K&T h j K =l 
3 Keujj 

= (d+ 1)^^ Ajj 

3 

<(d+ l)||tt|||maxAjj. 

3 

On the other hand, using the canonical basis vectors ej we have 

Amax(^) — e j Ae.j = Ajj, j = 1, . . . , N v i, 



and altogether we get (12). 

Using the same procedure for a diagonal scaling S = (sj) we obtain 



max( S - 2 Ajj) < A max (5-M5 _i ) < (d + 1) max(«TM«). 

3 3 



For the Jacobi preconditioning we have Sj = Ajj, which gives estimate (jl 3|) . □ 



Remark 4.2. Bound (13) can also be obtained by using the unassembled form of A as shown in |19| 
Sect. 3]. However, the analysis employed in Q15] cannot provide a lower bound on \ m \ n (S~ l AS~ l ) 
other than the trivial one, A m i ri (5'~ 1 J 4S'~ 1 ) > 0. 
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4.1 Geometric interpretation 



Although Lemma 4.1 gives a very tight bound on A max (^4), it does not provide any explanation on 
how the mesh or the diffusion matrix affect the conditioning. We now derive a bound on A max (A) 
in terms of mesh quantities and the diffusion matrix. 

Let Fx '■ K — > K be the affine mapping from the reference element K to the mesh element K, F' K 
the Jacobian matrix of Fk, jx the local index of 4>j on K and cj>j K = Fx ° 4>j K the corresponding 
basis function on K. 

Using the chain rule, we have 

Ajj = ^2 \K\ V<fij • DjfV^j 

= \K\{(F , K )- T V$ jK )-B K {(F' K )- T V$ jK ) 



< \ K \ (^r 1 ^^ 



/ \-T 



K<=u 



<CxY. \ R \ {F' k )- 1 ^k{F' k ) 



(14) 



where Ci = maxi jr= i r .. ) d4-i|| ^4>i K |||- Combining this with Lemma 4.1 yields 



A max (A) < (d+l)C.max ^ \ K \ 



(15) 



Remark 4.3. If we denote the maximum number of elements meeting at a mesh point by p max , then 
bound (15) implies 



Amax(-<4-) < 



■I \-T 



which is comparable to the estimates mostly found in the literature (e.g., [3 EJUS]). Note that both 
this bound and (15) are less tight than the bound (12). 



Remark 4.4. Bound (15) implies that the scaling 



s) = £ |JT| {F'^KiF^ 



l,...,N vi 



will also lead to bounds similar to (13) and those in Sect. |5| In general, this scaling is greater than 
the Jacobi preconditioning (cf. (14)) although they are equal in ID or for a mesh that is uniform in 
the metric specified by D _1 (cf. Sect. 4.1.3). 



4.1.1 Special case 



I 



For the simplest case of 



/, bound (15) has a rather simple interpretation. The quantity 



IK-F^-) ||2 can be bounded by the reciprocal of the in-diameter h m - mj K of K f!51 Lemma 5.1.2]. If 

- - ..I 1 1 

we denote the average aspect of K by Hk (i-e., Hk = \K\ d ), then we can rewrite (15) as 

2 / - \ 2 



A m ax(^4) 5; Ci max ^2 \ K 



1 



h 



min,_R" 



Ci max \K 



rf — 12 
d 



h K 



J 



h 



min,X 



The ratio t — is a measure of the aspect ratio of K. Thus, for the case of D = I, the largest 
eigenvalue of A is bounded by the maximum volume- weighted element aspect ratio of the mesh. This 
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is consistent with the observation by Shewchuk in [TH] where a detailed discussion on the relation 
between the largest eigenvalue of the stiffness matrix and the element aspect ratio is available for 
the case of D = / in d = 2 and d = 3 dimensions. 



For the general case D / /, on the other hand, it is more convenient to interpret bound (15) in 
terms of the mesh quality measures introduced in [TT]. We now proceed with this. 

4.1.2 Mesh quality measures 

The first measure is the alignment quality measure, which can be simply viewed as an equivalent to 
the aspect ratio of K in the metric specified by D^ 1 . It is defined as 



±tr((F' K )-^ K (F' K ) 
det((F^B K (F' K )-TU 



and measures how closely the principal directions of the circumscribed ellipsoid of K are aligned 
with the eigenvectors of Hk and the semi-lengths of the principal axes are proportional to the 
eigenvalues |15j . Notice that 

1 < Qah,o-i( K ) < °°. 

In particular, Qaiifi- 1 ^) = 1 implies that K is equilateral in the metric D^ 1 . 

The second measure is the equidistribution quality measure defined as the ratio of the average 
element volume to the volume of K, both measured in the metric specified by U^- 1 , 



Q eq p-i(K) = T gp- 1 (l(i) 



where | A^-i = \K\ det(D^) 2 is the volume of K with respect to D^ 1 and 

= £ \ K W ■ (17) 

The equidistribution quality measure satisfies 

< Q eq ,n-<K) < 00 and ^ £ QeqH-i( K ) = 1. 



Notice that 



a h = \K\det(B K )-% -»■ f det(D(oj)) » dx = |0| 



K&T h 

as the mesh is being refined. As a consequence, an can be considered as a constant. 
4.1.3 Geometric interpretation (general case) 

Using the quality measures we can rewrite the key factor ||(i 71 ^.) _1 Dx(i ? ft-)~ T ||2 as 
{F' K )- 1 ® K {F' K )- T n < tr UF' K )- l I} K {F' K ) 



2 

= dQ^r i (K)(\K\det(n K 



'all,] 



_> 



,.'N\d 
a 



2 2 

d 



Qil.iJn,-i( K )QeqM-i(K) 
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and therefore 



A max (A) <C(—Y max £ \K\ Q^. 1 B _ 1 (K)Q eq>0 - 1 (K) 



(18) 



Thus, A max (^4) is bounded by the maximum volume- weighted, combined alignment and equidistri- 
bution measure of the mesh in the metric D^- 1 . 

When a mesh is adapted to the coefficients of the BVP, i.e., it is uniform in the metric I 
will have the properties 



it 



and 



Qaup-i(K) = l, Q eq>0 -i(K) = l, VKeT h 



2 



< 



2 



(19) 
(20) 



Moreover, bound (18) will reduce to 



A ma x(vl) < CNd ItOn 



5 Smallest eigenvalue and condition number of the stiffness matrix 

The approach employed in this section was originally developed by Bank and Scott P for isotropic 
meshes. We generalize it here to arbitrary anisotropic meshes. 

Hereafter, we will use C as a generic constant which can have different values at different 
appearances but is independent of the mesh, the number of mesh elements, and the solution of the 
BVP. 

We start with bounds on A mm (^4). 

Lemma 5.1 (Smallest eigenvalue). The smallest eigenvalue of the stiffness matrix for the linear 
finite element approximation of BVP ([I]) is bounded from below by 



n, 



Amin(A) > Cd miQ N- 



1 +ln 



\K\ 



| if mini 

* (H) 



KeT h 



for d = 1 , 
ford = 2, 

for d > 3, 



(21) 



where \K\ = denotes the average element size. 

The smallest eigenvalue of the diagonally (Jacobi) preconditioned stiffness matrix is bounded from 
below by 



X^iS-'AS- 1 ) > CN~ 



1 



Nd ri 



E : 



\K\ 
~K\ 



-i 



for d = 1 



(22) 



and 



XminiS^AS- 1 ) > CN' 



\K\ {F' K )^ K {F' h 



K) 



X < 



1 + 



max \\(F')- l O K (F') 



In 



£ i^ilK^)- 

KeT h 



l \-T\\ 

k> 1 1 -a 



-1 



for d = 2, 
for d > 3. 



(23) 
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Proof. Since Sobolev's inequality is different for d = 1, d = 2 and d > 3 dimensions [HI Theorem 7.10], 
we treat these cases separately. 

Case d = 1: Let Cs be the constant associated with Sobolev's inequality. Using the inequality 
([2]), Sobolev's inequality, and the equivalence of the vector norms, 

u T Au > d m i n 

> d min C s l^l -1 sup \u h \ 2 
n 



= d m i n Cs 1 max u 2 
j J 

> d^nCsl^r 1 N- 1 \\ U \\l 



Therefore, A min (^) > Cd^N' 1 . 
With scaling, 



T ■ s 2 s~ 2 u 2 

v^S^AS^u > Cd min maxs7 2 u 2 > Cd mm 3 J \ 1 
In ID, 'V(pj K \k = an d therefore 



mm "«>a o j u,j ^_ <_vu. mm 2 — C d^ain 2' (^^0 



Using this in (24) gives (22). 

Case d = 2: Consider a set of not-all-zero non-negative numbers {ok, K G 7h} (to be determined 
later) and a finite number q > 2. Let Cp, Cg, and be the constants associated with Poincare's 
inequality, Sobolev's inequality, and the norm equivalence on K, respectively. Using Poincare's, 
Sobolev's and Holder's inequalities and the norm equivalence for Uh, we have 



u 1 Au = / Vuh ■ OVuh dx > d m i n \ Vuh\ H i 
Jn 



(«) 



> 



> 



dmin Cp 

l + Cp 



\ U h\\m(Q.) 



dminCpCs 1 M 2 



1 + Cp q 
dminCpCs 1 



\ u h\\Li(n) 



1+r/ , - U E KIIL W 



q-2 , g-2 2 



3 / „ \ 1 I \ q 

\1 



1 \ I 1 



dmmCpCs 1 | ;, j 1 | ; .' 1 ( V 



dminCpCs 1 
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9-2 
9 



1 + Cp g \ ReTh / XeTh 



g-2 
1 



dminCpCsl I \ - \ ir-i=ii- 1 

1 + op g ^ eTh y A , eTh 

. dminCpCsC^ 1 / ^ ^ 17^1-n ||2 
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dminCpCsCfi 1 



1 + Cp q 



E< 



q-2 
<3 



9-2 

if 



E n j E a A - • 



The choice ax = \K\ " yields 



u T Au > C(L 



E l*l 



'q-2 



q-2 
i 



and therefore 



A min (^) > Cd^q- 1 [ ^ |K| 
yK&T h 



q-2 



2 

"?-2 



q-2 

<3 



> CtW* (iV \K n 

= Cd min N- 1 Iq-HNlK^yq 



q-2 



(25) 



The largest lower bound on (25) is obtained for q = max {2, |ln (N \K m i n \)\) with 

C 



-l 



(iV|^ min |H > 



l + |ln(iV|iC min |)| 



The choice q = 2 is viewed as the limiting case as q — > 2 + . Estimate (21) follows from this, (25) 
and the definition of the average element size. 
With scaling, we have 



vfS^AS^u > Cd n 



1 



q 

E a K 2 



q-2 
q 



E""I s j 2 E q k\k 



J 3 

j Keu 



For the Jacobi preconditioning sj = Ajj = Y^Keu>j ' ^K^^j we choose 

-2 d+1 -2 d+1 

a K = \K\~r V<f> iK -B K Vct> iK = \K\-T £ V$ iK ■ {F' K )- l B K {F' K )- T V4> iK , 

i K =l 



i K =l 



which gives 



and 



x „ 2 

> 0!R \K\ i > 1 



q- 2 



olk < {d + l)C,\K\— (F- 



" k )- 1 ^k{F k 



■i \-T 



where = maxj if= i v . i) d_|_i||V0j K || 2 . With these and choosing the value for the index q in a similar 
manner as for the case without scaling we obtain (23). 

Case d > 3: This case is very similar to case d = 2. Again, from (|2j), Poincare's, Sobolev's and 
Holder's inequalities and the norm equivalence for Uh, we have 



U T Au = / VUh ■ OVUh dx > dminl Vu/jjjyl 

Jn 

^ dminCp || || 2 

- Y+c~p~ ™ H1 w 



(fi) 
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. dminCpCs n n 2 

> : , „ \\u h \\ 2d_ 

1 + Cp Ld-2(Q) 
dminCpCs 



I i f< I Z^, W U h\i 

1 + Cp \K&Th L d ~HK), 



2d 
| d-2 
2d 



d-2 
d 



2 2 d-2 

d \ d ( d \ d ( 2d_ ^ d 



E <4 E E Kir. 



2 

2d 



. dminCpCs [ I I n i|2 

> — i -p, / "ft- > OLK \\Uh\\ 2d 

1 + Cp \ *—L I *—t Ld=2(K) 

1 ^^P \K<=T h J K&T h L (K) 



dminCpCs I v-^ § \ \^ |r ,,%2 „ 

: , „ Z_v a K\ OLK\K\ * \\U h \\ 2d_ „ 

1 + Cp Uer, Ke% Ld ~ 2 ^ 



d 



1 + Cp 

dmmCpCsCfr 
1 + Cp 



_ d-2 

The choice ax = \K\ d gives 



d min C P CsCft I x - , | \ - i/'i : '7 J ii i|2 



K KeT h I KeT h 



E 4 E n i E a K\K\ * 
,K&T h I j Keuj 



u T Au>Cd min [ E ] E 



u). 



Estimate (21 ) follows from this and the definition of the average element size. 
The bound for the scaled stiffness matrix is obtained by choosing 

2 d+1 

a K = \K\S V0 iK • (F^)- 1 D^(F^)" T V^ K . □ 



, 2 

i K =l 



Combining Lemma 4.1 estimate (15) and Lemma 5.1 we obtain upper bounds on the condition 
number of the stiffness matrix and the scaled stiffness matrix. 

Theorem 5.2 (Condition number of the stiffness matrix). The condition number of the stiffness 
matrix for the linear finite element approximation of BVP Q is bounded by 



1 \K\ 
k(A)<CN 2 - — max £^ BjfJ— [, for d = 1 (26) 



and 



(l— - 
— -max Y, \ R \ ( F k)~ l ®K(Fk)~ T 



1+lnpg^, ford = 2, 

, - d-2\ 1 (27) 
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The condition number of the diagonally (Jacobi) preconditioned stiffness matrix is bounded by 



uis-'as- 1 ) < Civ 2 ^— J2 ° K wr f° rd = 1 



and 



niS^AS- 1 ) < CN< 



• V ' /|!,il; A T T h 



£ \K\ {F ! K r^ K {F' K ) 



I \-T 




1 + 



111 



max \\{F' K )-^ K {F' K )- T \ 



\k\\\(f' k )-ih k (f> k )-t\ 

xer h 



(28) 



, ford = 2, 



(29) 



for d > 3. 



5.1 Geometric interpretation 

We now study the geometric interpretation of the bounds for the condition number. 



5.1.1 Without scaling 



Bounds (26) and (27) contain three factors, a base bound CNd } a factor reflecting the effects of 
the mesh nonuniformity measured in the metric D" 1 (mesh 1$ -nonuniformity) , and, if d > 2, a 
factor reflecting the effects of the mesh nonuniformity in volume measured in the Euclidean metric 
(volume-nonuniformity). 

2 

The first factor N a corresponds to the condition number of the stiffness matrix for the Laplacian 



operator on a uniform mesh (cf. Special Case 5.3 below). 
The second factor 

{F' K y^ K (F' K ) 



max 

j 



K 



reflects the effects of the mesh D-nonuniformity and can be understood as a volume-weighted, 
combined alignment and equidistribution quality measure of the mesh with respect to D _1 (cf. 



Sect. 4.1.3). 



The third factor in (27) is 



1 +ln 



\K\ 



|-Kmin| ' 

N rhr. W 



for d = 2, 
for d > 3. 



KeT h 



It measures the effects of the mesh volume-nonuniformity (measured in the Euclidean metric) on 
the condition number. Notice that there is no effect in ID and in 2D it is minimal. In d > 3 
dimensions the factor is proportional to the average of \K\ d over all elements. This is a significant 
improvement in comparison with previously available estimates which are proportional to |-fiTmin| — 



HI- 



5.1.2 With scaling 



Bounds (28) and (29) for the scaled stiffness matrix have the same base bound as without scaling. 
Hence, diagonal scaling has no effect on the condition number when the mesh is uniform and D = I. 
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Unlike (27), bounds (28) and (29) do not have the third factor which involves only the element 
volume (in comparison to the second factor which couples (F K )~ l with D~ x ). In this sense, a 
properly chosen diagonal scaling can eliminate the effects of the mesh volume-nonuniformity on 
the condition number. Moreover, scaling can also significantly reduce the effects of the mesh D- 
nonuniformity. Indeed, the factors in (28) and (29) that couple {F K ) with D _1 are asymptotically 
the L max{1 'i } (ft) norm of \\(F' K y 1 B K (F' K )- T \\ 2 whereas the corresponding factors in (|26j) and ((27]) 
are basically the maximum norm. 

Furthermore, the D-related factor in (29) for d > 2 can be rewritten in terms of the alignment 
quality measure Qaim- 1 from Sect. 4.1.2|as 



E \K\ (F' K )~ l H) k{F'k)~ T l<dt £ Q^(K) dei(B K )l 



(30) 



KeT h 



Thus, the dependence of this D-related factor on the element volume is also mild: both Q a n jj-i (K) 
and Bk (the average of D over K) are invariant under the scaling transformation of K. 

The following special cases are instructional to understand the interplay of the factors for different 
types of meshes. 



Special Case 5.3 (Uniform meshes). For a uniform mesh and D = I, bounds (26)-(29) yield 

k(A)<CnI and ^S^AS' 1 ) < Cm, 



which is the base bound. Hence, the diagonal scaling has no effect on the condition number when 
the mesh is uniform and D = /. 



Special Case 5.4 (Isotropic meshes, 



/, d > 2). For an isotropic mesh and D = I, 



\K\ 



h K and 



'F\ 



K) 



-r 



K 



Therefore, 



1 

N 



J2 \K\ {F' K r^ K {F' h 



-T 



KJ 



KeT h 



1 < _ 

2 ~ N 



E f4ch 



d u-d 
K 



Ken 



and bound ( 29 ) reduces to 

Kis^AS- 1 ) < cm 



l + ln 



1^1 



for d = 2, 
for d > 3, 



(31) 



which is precisely the result of Bank and Scott [2, Theorems 4.2 and 5.2]. In this case, the diagonal 
scaling becomes 



(A, 



03 



E 4~ 2 

\Keujj / 



where hj denotes the average length of the elements around the j vertex. This scaling is equivalent 
to the change of basis functions 



2-d 
2 



<f>i -»■ h j 2 <j>j 



which is used in [2J Example 2.1] 
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tl-1 



Special Case 5.5 (Uniform meshes with respect to 
to D _1 , i.e., coefficient adaptive, we have properties (ll 91) and (20). Bounds (26)-(29) reduce to 



). For a mesh that is uniform with respect 



k{A) < 



CCN>max|) 



C 



dr, 



N 
o-h 

N 



1, 

1 +ln 



1*1 

| if min | ' 



for d 
for d 



V Mr. KM 



N KeT h 
for d> 1, 



for d > 3, 



where cr^ is defined in (17) and corresponds to the volume of the domain in the metric specified by 
D _1 . Thus, the condition number of the scaled stiffness matrix for a coefficient adaptive mesh has 
the optimal order of 0{Nd). 

Special Case 5.6 (Aligned meshes, d > 2). For meshes aligned with the diffusion matrix but not 
necessarily fully coefficient adaptive (i.e., isotropic but not uniform with respect to D _1 ) we have 



Q. 



all,} 



-x(K) = l but Q eq n-i(K)^l. 



From (30), bound (29) becomes 



KiS^AS- 1 ) < C 



2 

Nd 

^min 



N 



E det ( 



KeT h 




1 + ln 



1*1 

l^minl ' 



for d = 2, 
for d > 3. 



Aside from the term depending on det(D), this bound is equivalent to (31). Hence, the diagonal 
scaling almost eliminates the effects of the mesh on the condition number for B-aligned meshes. 

Special Case 5.7 (General M-uniform meshes). Finally, let us consider general M-uniform meshes, 
i.e., meshes that are uniform in the metric specified by a given metric tensor M which does not 
necessarily correspond to D _1 . In the context of mesh adaptation, an adaptive mesh is typically 
generated based on some estimate of the solution error and the associated metric tensor M is 
solution dependent. Thus, it is of interest to know what the impact of a given M on the conditioning 
of the stiffness matrix is. Recall [13^ that an M-uniform mesh satisfies 



{f'kT t {f' k ) 



N 

0~h,M 



M K , 



where Mk is some average of M on K and ah,M is defined as in (17) but with 
We have 



replaced by M 1 . 



K ) 



{F'k) 



and therefore 



k^AS- 1 ) < 



C 



-T 



N 



< 



N 

0~h,M 



\\M K B 



K\\ 2 



E \K\\\M K B K 




Hence, the bound on the condition number after diagonal scaling for an M-uniform mesh depends 
only on the volume- weighted average of HMxOif or, asymptotically, the L d l 2 norm of ||MD|L. 
For many problems such as those having boundary layers and shock waves, mesh elements are 
typically concentrated in a small portion of the physical domain. In that situation, we would expect 
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that M differs significantly from D only in small regions. As a consequence, the volume-weighted 
average of H-Mr-B^H^ over the whole domain may remain small and therefore the condition number 
of the scaled stiffness matrix for anisotropic adaptive meshes does not necessarily increase as much 
as generally feared. 



This effect can be observed in Examples 6.2 and 6.3. Figures 3a and 4a show that the effects 



of anisotropic adaptation are completely neutralized by the diagonal scaling when the number of 
anisotropic elements is small in comparison to N. 



6 Numerical experiments 



In this section we present numerical results for a selection of one-, two-, and three-dimensional 
examples to illustrate our theoretical findings. 

Note that all bounds on the smallest eigenvalue contain a constant C. We obtain its value by 
calibrating the bound for A m i n (5 _1 74S' _1 ) with Delaunay (Example 6.4) or uniform meshes (all 
other examples) through comparing the exact and estimated values. For the largest eigenvalue we 



use explicit bounds (12) and (13). 



First, we give examples with predefined meshes to demonstrate the influence of the number 
and shape of mesh elements on the condition number of the stiffness matrix and to verify the 
improvement achieved with the diagonal scaling. For the tests, we employ the Laplace operator (i.e. 
B = I) and a mesh on the unit interval, square, and cube, for ID, 2D, and 3D, respectively. 

Example 6.1 {d = 1, D = /, Chebyshev nodes). For a simple one-dimensional example we choose 
a mesh given by Chebyshev nodes in the interval [0, 1], 



1 



cos 




,N — 1. 



(32) 



The exact condition number of the stiffness matrix and its estimates (26) and (28) are shown in 



Figs, [la] (without scaling) and |lb| (with scaling) while those for the extreme eigenvalues and their 



estimates are given in Figs, lc (without scaling) and |ld| (with scaling). 



Figure [la] shows that the estimate (26) is much sharper than the standard estimate with A m i n (A) oc 
|-^min|- The former has the same asymptotic order as the exact value as N increases, whereas the 
latter is too pessimistic and has a higher asymptotic order. The difference is caused by the estimate 



of the smallest eigenvalue (Fig. lc). Notice that the estimates on the largest eigenvalue are very 



tight, both for the scaled and the unsealed cases. 

The results clearly show the benefits of diagonal scaling: the order for the condition number of 
the scaled stiffness matrix in Fig. lb is 0(N 2 In N), which is almost the same as for uniform meshes, 



whereas that without scaling in Fig. [lalis 0(iV 3 ). It can be shown analytically that the orders of 



the nonuniformity factors in (26) and ( |28[ ) for the Chebyshev nodes defined with (32) are O(N) and 
C(ln N) and those of the corresponding condition numbers are 0(N 3 ) and 0{N 2 In N). 

Thus, the numerical and theoretical results are consistent and the improvement by diagonal 
scaling from the maximum norm to the L 2 norm is significant in this example. 



/, anisotropic elements in a unit square). For this 2D example we use a 

First, we 



2a 



Example 6.2 (d = 2, ] 

mesh for the unit square [0, 1] x [0, 1] with 0(7V 1//2 ) skew elements, as shown in Fig. 
fix the maximum aspect ratio at 125 : 1 and increase N to verify the dependence of the condition 
number on N (Fig. 3a). Then, we fix N at 20,000 and change the maximum aspect ratio of the 



mesh elements to investigate the dependence of the conditioning on the mesh shape (Fig. 3b ) 



Figure 3a shows the averaging effect of the diagonal scaling: the scaling significantly reduces 
the condition number and, when N becomes large enough, the conditioning of a scaled system is 
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Figure 1: Example 6.1 Exact and estimated condition number and eigenvalues of the stiffness 
matrix as a function of N (d= 1). 



comparable to the condition number on a uniform mesh. Moreover, the estimated value of the 
condition number with or without scaling has the same order as the exact value as N increases. 



Figure 3b provides a good numerical validation of (27), namely that the condition number of 
the unsealed stiffness matrix is linearly proportional to the largest aspect raticj^] With scaling, the 
condition number is still increasing with an increasing aspect ratio, since the average aspect ratio 



(in accordance to (29)) is also increasing. Nevertheless, the condition number after scaling is smaller 
by a factor of 10. 



Figure 3b also shows that our estimate of the condition number with scaling has the same (linear) 
order as the exact value as the maximum aspect ratio increases, whereas the bounds for the unsealed 
case has a slightly higher order. This indicates that the estimation can be further improved. 



As for the estimates on the extreme eigenvalues, the results are mainly the same as in Example 6.1 
For this reason, we omit them in 2D and 3D to save space. 



Example 6.3 (d = 3; anisotropic elements in a unit cube). In this example, we repeat the same 
test setting as in Example 6.2 fixed anisotropy (25 : 25 : 1) with increasing number of elements 
(Fig. 4a) and a fixed N = 29,478 paired with the changing anisotropy of the mesh (Fig. 4b). The 
results shown in Fig. [4] are essentially the same as in 2D. Since the mesh used in this example has 



In 2D with D = I, the nonuniformity term in (271 is equivalent to the aspect ratio. 
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(a) ©(TV 1 / 2 ) skew elements. 



(b) 0(N 2 / 3 ) skew elements. 



Figure 2: Predefined meshes for (kfo Example 6.2 and <von Example 
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Figure 3: Example 6.2 Condition number before and after scaling for a predefined 2D mesh (Fig. 2a) 
as a function of ([a]) the number of mesh elements and (|b]) the maximum element aspect 
ratio. 
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Condition number before and after scaling for a predefined 3D mesh (Fig. 2b ) 
as a function of ([a]) the number of mesh elements and (|b]) the maximum element aspect 
ratio. 
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(b) Numerical solution. 



Figure 5: Example 6.4: Boundary conditions and a numerical solution. 



a larger share of skew elements {0{N^ 1 ^)) than the mesh used in Example 6.2 (0(iV -1 / 2 )), it is 
reasonable to expect that the averaging effect of diagonal scaling is less effective. This can be seen 
in Fig. [4] where the exact condition numbers with and without scaling stay closer than in Fig. [3} 

Figure [4] shows that the bounds on the condition number with and without scaling have the 
same asymptotic order as the exact values as N increases. However, they have higher orders as the 
maximum aspect ratio increases for a fixed N . As in the previous example, this indicates that the 
estimation can be further improved. 

In the next example, we consider an adaptive finite element solution of an anisotropic diffusion 
problem with different meshes. 

Example 6.4 {d = 2, adaptive anisotropic meshes). Consider an anisotropic diffusion problem 
studied in |16| [T7] . It takes the form of BVP ([TJ but with a non-homogeneous Dirichlet boundary 
condition. The domain and its outer and inner boundaries T ou t and Ti n are shown in Fig. [5aJ The 
coefficients of the BVP, the right-hand side and the boundary data are given by 



cos ip — sin ip 
sin ip cos ip 



f = 0mn = [0,l} 2 \ 



1000 
1 

4 5^ 2 ' 



cos ip sin ip 
— sin ip cos ip 



g = on r 



outi 



9 



ip = 7r sin x cosy, 



2 on r,- 



(33) 



We employ an adaptive finite element algorithm from |14l IT6] to compute the numerical solution 
and adaptive meshes. The algorithm utilizes the M- uniform mesh approach, i.e., meshes are 
generated as quasi-uniform in a given metric M. For the mesh generation we use the bidimensional 
anisotropic mesh generator [2]. 

A Delaunay mesh — our first example (Fig. 6a) — is M- uniform (or M-quasi- uniform) with respect 
to M = I. The second mesh (Fig. |6b[ ) is purely coefficient-adaptive and is defined as an M-uniform 
mesh with respect to D , i.e., M = D . The third mesh (Fig. 6c) is a purely solution-adaptive 
mesh where M = M{uh) depends on the numerical solution Uh (or, more precisely, on the hierarchical 



basis error estimate e^). The fourth mesh (Fig. 6d) represents a combination of adaptation to both 
the solution and the coefficients of the problem and the metric is defined as M = ^(e/^B -1 , where 
9 {eh) is a scalar function depending on the error estimator e\ t . With such choice the shape of mesh 
elements is determined by the diffusion matrix while the size is controlled by the estimate of the 
solution error. 

From Fig. [6] we can see that the smallest condition number among all four meshes is with the purely 
coefficient-adaptive mesh (Fig. 6b), which is consistent with Special Case 5.5 The conditioning is 
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Figure 6: Example 6.4 condition number in dependence of N. 
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better than in the case of a quasi-uniform mesh (Fig. 6a), confirming the observation that, depending 
on the problem, a quasi-uniform mesh is not necessarily the best mesh from the conditioning point 
of view. For both cases, diagonal scaling does not improve the condition number significantly. This 
is expected since both meshes are almost volume-uniform. To explain why the mesh in Fig. 6b is 



(almost) volume-uniform, we recall from (16) and (19) that an M-uniform mesh with respect to 
satisfies 

VK e r h . 



\K\ 



det(D 



'K) 

1000. Thus, 



\K\ = const. 



The diffusion matrix D in (33) satisfies det(B 

The largest condition number is in the case of the purely solution-adaptive mesh (Fig. 6c ) . This 
is because the mesh is not volume- uniform and its elements are not aligned with D _1 . Since the 
mesh is far from being uniform in size, scaling will have a significant impact, as it can be verified 
in Fig. 6c the condition number after the scaling is even smaller than the condition number with 
Delauney meshes. 

Conditioning with a mesh that is both coefficient- and solution-adaptive (Fig. 6d ) is not as good 
as in the case of the purely coefficient-adaptive mesh but better than in the case of the purely 
adaptive and Delaunay meshes. 

In all four cases we observe that the developed estimates for the condition number of the stiffness 
matrix are reasonably tight and have the same order as the exact values as N increases for both 
unsealed and scaled cases. 



7 Summary and conclusions 

Mass matrix 

Our new estimate ^ of the condition number of the Galerkin mass matrix is tight within a factor of 
{d+ 2) from both above and below for any mesh with no assumptions on mesh regularity or topology. 
It this sense, it is optimal and truly anisotropic. 

Stiffness matrix 



Lemma |4.1| provides an estimate of the largest eigenvalue of the stiffness matrix which is simple to 
compute and is tight within a factor of (d+ 1) from both above and below for any mesh. This is in 
contrast to many existing estimates which are proportional to the maximal number of elements 
meeting at a mesh point. 



New bounds (21 )-(23) on the smallest eigenvalue and (26)-(29) on the condition number of the 
stiffness matrix are a significant improvement in comparison to the previously available estimates. 

First, the new bounds show that the conditioning of the stiffness matrix with an arbitrary 
(anisotropic) mesh is much better than generally assumed, especially for d = 1 and d = 2. 

Second, the new bounds are truly anisotropic and valid for any mesh since no assumptions on the 
mesh regularity were made. 



Third, bounds (26) and (27) reveal what affects the conditioning. There are three factors. The 

2 

first (base) factor CNd describes the direct dependence of the condition number on the number 
of mesh elements and corresponds to the condition number for the Laplace operator on a uniform 
mesh. The second factor describes the effects of the mesh D _1 -nonuniformity, i.e., the interplay 
between the shape and size of mesh elements and the coefficients of the BVP. It is 0(1) for a 



coefficient-adaptive mesh, i.e., a mesh satisfying (19). The third factor measures how the mesh 
volume-nonuniformity further affects the condition number. It has no effect in ID, a minimal one 
in 2D, and a substantial effect in 3D and higher dimensions. This means that even if the mesh is 
coefficient-adaptive and the second factor is 0(1), the mesh volume-nonuniformity can still have a 
significant impact on the condition number for d > 3. 
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Fourth, a simple diagonal scaling, such as the Jacobi preconditioning, can significantly improve 



the conditioning. Bound (29) for the condition number after scaling does not contain the factor 
for the mesh volume-nonuniformity. As a consequence, for a coefficient-adaptive mesh, this bound 
reduces to the base factor CNd. In this sense, diagonal scaling eliminates the effects of the mesh 
volume-nonuniformity. It can also significantly reduce the effects of the mesh nonuniformity with 
respect to D _1 : the influence reduces essentially from the maximum norm to the L max l 1 >2} norm of 
\\(F' K )-^ K (F' K )- T \\ 2 . 

Moreover, for a preconditioner that is invariant to diagonal scaling it follows that the condition 
number of the preconditioned stiffness matrix is typically smaller than k(S~ 1 AS~ 1 ) which in turn 



has a much lower bound than k(A) (cf. ( |27[ ) and (29)). For example, consider an incomplete 
Cholesky decomposition of A, 

A = LL T + E. 

It follows that 

S^AS- 1 = ( S- l L) ( S^lY + S^ES- 1 



is actually an incomplete Cholesky decomposition of S 1 AS 1 since S 1 L has the same sparsity 
pattern as L. Then from the identity 



Lr x ALr T = (s^l) 1 (s-'as- 1 ) (s-'l) 



we see that the preconditioned matrix of A with preconditioner L is equivalent to the preconditioned 
matrix of S~ 1 AS~ 1 with preconditioner S~ 1 L. As a result, the performance of the preconditioning 
technique on A is the same as that on S~ 1 AS~ 1 which has a much smaller condition number than A. 
Although there is no estimate yet on the condition number of the preconditioned system, the above 
observation may provide a partial explanation for the good performance of ILU preconditioners 
with anisotropic meshes observed in 



Numerical experiments (Figs. 3a and 4a) indicate that although the new bounds have the same 
order as the exact value as the number of elements increases, they may have higher asymptotic orders 
than the exact value as the element aspect ratio increases. These may deserve further investigations. 

Finally, we would like to point out that although the study in this paper has been done specifically 
for the linear finite element discretization, the approach can be generalized for higher order finite 
elements without major modifications. 
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